Long-time self-diffusion of Brownian Gaussian-core particles 
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Using extensive Brownian dynamics computer simulations, the long-time self-diffusion coefficient 
is calculated for Gaussian-core particles as a function of the number density. Both spherical and rod- 
like particles interacting via Gaussian segments are considered. For increasing concentration we find 
that the translational self-diffusion behaves non-monotonically reflecting the structural reentrance 
effect in the equilibrium phase diagram. Both in the limits of zero and infinite concentration, it 
approaches its short-time value. The microscopic Medina-Noyola theory qualitatively accounts for 
the translational long-time diffusion. The long-time orientational diffusion coefficient for Gaussian 
rods, on the other hand, remains very close to its short-time counterpart for any density. Some 
implications of the weak translation-rotation coupling for ultrasoft rods are discussed. 

PACS numbers: 66.10.Cb; 61.20.Ja; 82.70.Dd 



I. INTRODUCTION 

Particles interacting via penetrable pair potentials ex- 
hibit fascinating new clustering and reentrance effects 
[3, 0, H, 0, 0] which are absent for diverging potentials 
such as hard spheres and inverse-power potentials. A 
well-studied model for a penetrable interaction is a Gaus- 
sian potential [!, 0, 0, 0> 01 which mimics the effective 
interactions between two polymer coils in a good solvent 
[l(| and applies also to dendrimer solutions [ll], Qjl, [l3| . 
This potential can be generalized towards a Gaussian- 
segment model for rod-like particles in order to describe 
bottlebrush polymers with a stiff backbone [l4[ , see also 
[l5l |. In the Gaussian-core system, two particles pay a 
finite energy penalty if they are sitting on top of each 
other. If they overlap completely there is no repulsive 
force any longer. 

In the present paper we focus on equilibrium dynamical 
correlations of Gaussian Brownian fluids. In particular 
the long-time self-diffusion coefficient is simulated as a 
function of the particle density. Recently the dynami- 
cal behavior of spherical Gaussian particles has been ex- 
plored by molecular dynamics studies [l6j which is suit- 
able for polymer melts but neglects the hydrodynamic 
friction of a solvent. Here, we consider solutions of col- 
loidal or polymeric particles and therefore overdamped 
Brownian dynamics is appropriate where the friction of 
the solvent is included. We further study the long-time 
translational and orientational self-diffusion of Gaussian 
segment rods in the isotropic phase as a function of rod 
concentration. 

As a result, we find that the long-time self-diffusion 
coefficient behaves non-monotonically with density, both 
for spheres and rods. For zero density (i.e. the sin- 
gle particle limit) the long-time self-diffusion coefficient 
is clearly identical to the short-time diffusion constant 
which is entirely dominated by solvent friction. What is 
less obvious is that for high densities with multiple over- 
lap of particles the self-diffusion again tends to its short- 
time counterpart. In fact in the limit of very high densi- 



ties, a Gaussian particle feels many neighbors around a 
distance where the Gaussian potential has its inflection 
point and these give rise to a diverging number of inter- 
action kicks. Therefore one could have expected a higher 
diffusion coefficient than the short-time value. However, 
we show that correlations between the neighboring par- 
ticles enforce a normal diffusive behavior of an effective 
ideal gas in this limit. 

Between these two extreme limits, for finite densities, 
the long-time self-diffusion coefficient is smaller than its 
short-time counterpart. The minimal value is roughly 
at the point of maximal fluid structure. A similar non- 
monotonic behavior has been found for molecular dynam- 
ics [16| where the ballistic limit of zero-density leads, of 
course, to a diverging long-time self-diffusion coefficient. 
Furthermore we find that the long-time orientational self- 
diffusion coefficient practically coincides with its short- 
time behavior. The latter fact implies that there is no 
significant translation-rotation coupling in the Gaussian 
segment model for penetrable rods. 

We compare our simulation data with the microscopic 
theory of Medina-Noyola [l7], EH which relates the long- 
time self-diffusion coefficient to the fluid pair correlation 
and find qualitative agreement for spheres. The same 
holds for the translational diffusion of rods if the theory of 
Medina-Noyola [lTj is applied to the translational degrees 
of freedom alone. 

The paper is organized as follows: in Section II we 
describe in detail the procedure of the Brownian dynam- 
ics computer simulation. Results for the long-time self 
diffusion are presented and discussed in Section III. Fi- 
nally Section IV is devoted to more general remarks and 
conclusions. 



II. BROWNIAN DYNAMICS COMPUTER 
SIMULATIONS 

The Brownian dynamics (BD) simulations are based 
on a finite-difference integration of the overdamped 
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Langevin equations for N interacting anisometric parti- 
cles in three dimensions [l9|, |2(J. The trajectory of each 
particle i is characterized by its position (t) and orien- 
tation Cji(t) at time t. If hydrodynamic interactions are 
neglected, the update equation for the position of particle 
i can be written in the following way 



Ti(t+At) =Ti(t)- 



At 



Dj-F 4 (t) + Ar ! + 0{(At) 2 }, (1) 




FIG. 1: Soft rod of length L composed of Ns = 3 segments 
with intersegment distance A. 



with ksT the thermal energy and F$(t) the total force 
acting on the center-of-mass. The latter is derived from 
the pair potential which will be specified later. Fur- 
thermore, Dq represents the short-time diffusion ten- 
sor which in case of uniaxially symmetric particles (e.g. 
cylinders) can be cast into the form 



(2) 



in terms of the translational diffusion coefficients parallel 
(Dq) and perpendicular (Dq) to the particle axis, with 
I the unit tensor and (g) a dyadic product. 

The contribution Ar^ denotes a random displace- 
ment of the particle due to collisions with the solvent 
molecules. Similar to Eq. (J2), it is convenient to decom- 
pose it into contributions parallel and perpendicular to 
the particle axis. Introducing two orthogonal unit vec- 
tors, Bu and eai, perpendicular to u>j we can express the 
noise term in Eq. ([T]) as 

Ai-i = Ar'lfa>i(t) + Ar^enO) + Arf 2) e 2i (t). (3) 

Here, Ar" and Ar^ 2 * represent Gaussian random dis- 
placements parallel and perpendicular to the symmetry 
axis. Both stochastic quantities have zero mean and their 
variance is 2D\At and 2-D^At, respectively. 

The orientational update equation for uii(t) reads 

tbi(t+£t) = u i (t)+-4^D$T i {t)xu i (t)+ALj i +0{(At) 2 }. 
k B l 

(4) 

Here, Dq denotes the short-time rotational diffusion co- 
efficient and Ti(t) the total center-of-mass torque acting 
on particle i. The noise contribution, 



AQi = xiiu(t) + x 2 e 2i (t), 



(5) 



is generated by means of two uncorrelated random Gaus- 
sian numbers, X\ and x 2 , both with zero mean and 
variance 2Dq At. After each step the new orientations 
(2>i(t + At) have to be renormalized to ensure that |cDj| = 1 
at all times. 

Obviously, for spherical particles the translational 
Brownian motion is completely decoupled from the ori- 
entations and the translational diffusion tensor Eq. @ 
becomes diagonal, i.e. = DqI. In this case, we need 
only consider the update equation for the positional co- 
ordinates Eq. ([T]). All update equations are exact up 



to order O(At) which suffices for the present purpose, 
provided At is chosen small enough. For a detailed dis- 
cussion of a second order update algorithm the reader is 
referred to Ref. [20| . 

The pair potential of the particles is given by an ultra- 
soft Gaussian potential. For spherically symmetric par- 
ticles we have: 



V2(r) = eexp [—(r/a)" 



(6) 



where a is the potential range which will henceforth 
serve as our unit of length. The amplitude is fixed at 
e = 5k B T. According to 3, the associated reduced tem- 
perature T* — ksT/e = 0.2 is much higher than the 
upper freezing temperature T* = 0.01 which guarantees 
a stable fluid state at any density. 

Apart from ultrasoft spheres we will also consider sys- 
tems of Gaussian rods |14| . Each rod has length L and is 
composed of Ns spherical segments placed at equidistant 
positions along the rod main axis with segment-segment 
distance A = L/(Ns — 1). The interaction potential be- 
tween two segments from different rods is again a Gaus- 
sian. The pair potential between two rods i and j is then 
given by a sum over all segment interactions: 



V2(Ti,Tj;U!i,Wj) 



K K 
a=-K fi=-K 



expHlr^l/a) 2 ], (7) 



where K = (Ns — l)/2 and r a p = (r< + aA&i) — (r^ + 
f3Aujj) the distance between segment a on rod i and f3 
on rod j (i j). We will consider slightly anisometric 
rods with Ag = 3 segments and L — 2a (see Fig. [1] for 
a sketch). Furthermore the segment- segment potential 
amplitude is e = bk B T . The short-time diffusion coeffi- 
cients of the rods depend on one-particle hydrodynamic 
effects. For these, we take the analytical results obtained 
for hard ellipsoids of length L and aspect-ratio p > 1 
reported by Tirado and co-workers pl| : 

D^ = — -2-Qnp + 0.839 + 0.185/p + 0.233/p 2 ), (8) 
2p 

D U = ^L(i np _ 0.207 + 0.980/p-0.133/p 2 ), (9) 
18 D T 

Dq = — — 2-(lnp - 0.662 + 0.917/p - 0.050/p 2 fl;0) 

with Dq = kgT '/&nrj s a the short-time diffusion con- 
stant of a sphere with radius a and n s the shear vis- 
cosity of the solvent. In the above, we have implicitly 
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identified L — pa, with p the hydrodynamic aspect-ratio 
of the Gaussian rods. In order to enforce a substantial 
translation-rotation coupling we take a value p = 5 which 
is larger than the interaction aspect-ratio L/a = 2. 

A natural unit of time is the Brownian time tb = 
a 2 1 Dq defined as the typical time a Gaussian particle 
needs to diffuse over a distance comparable to its own 
dimension. Let us further introduce Dq = 1 1 § dujT)J \ | , 
the isotropic orientational average of the diffusion tensor 
Eq. ©. For the spheres Dq = Dq while for the rods 



di = ;d« 



=:D" 



(11) 



With this result, we may compute the ratio of the single- 
rod rotational and translational relaxation times, i.e. 
Tq R /tq T ' = D%/p 2 a 2 D§ ==! 0.264, showing that the short- 
time orientational dynamics is much faster than the 
translational. The quantity Dq also provides the natural 
scale for the long-time translational self-diffusion coeffi- 
cient Dj, defined as 



rp 11 

Dj = lim 

L t^oo Qt N 



Mi)-r*(0)) J 



(12) 



where (• • • } is a canonical average. An alternative defi- 
nition is provided by a differential expression 



D 



lim 



Id 1 



t^oo 6 dt N 



>,(t)-r 4 (0)r 



(13) 



Both expressions should in principle yield identical re- 
sults in the long-time limit. In practice however, they 
will differ slightly and the difference can be used to as- 
sess the error in Dj ■ 

Following Ref. [2(| we can define the long-time rota- 
tional diffusion coefficient Df; as follows 



Dl 



lim 



W n (t) 



t 



(14) 



where W n (t) is an orientational correlation function mea- 
suring the mean-square displacement on the unit sphere. 
It is given by 



W n (t) 



1 



n(n + 1) 



In <P„(w(t).(2>(0))) 



(15) 



with V n a Legendre polynomial. Similar to Eq. (|13p . we 
may also employ the differential analogue. If the rota- 
tional motion is a diffusion process on the unit sphere, 
the dynamics is captured by the Debye diffusion equation 
which predicts to be independent of n [2(| ■ 

In our simulations we used a cubic simulation box of 
volume V with periodic boundary conditions in all three 
directions. The number of particles depends on the num- 
ber density. If we choose to cut off all segment- segment 
pair interaction for which v 2 /e < f™*/ 6 the length Lb 
of the simulation box must be at least twice the cor- 
responding cutoff range and so Lb /a > — 2 ln^"'/ 6 ]- 
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FIG. 2: Long-time translational self-diffusion coefficient di- 
vided by its low-density limit D \ / Dq versus the density pa 3 
(on a log scale) for Gaussian spheres from BD simulations 
and the theory of Medina-Noyola (MN). On the right vertical 
axis, filled dots give the amplitude of the first maximum 
of the equilibrium pair correlation function (/2(f)- 



For a given density the number of particles must there- 
fore obey N > (pa 3 )L 3 Bl while imposing a minimum of 
N = 500 at small densities. Using v™ % /e 



at pa — 10 requires N = 4000 and one at pa 



25 



N = 13000. The time step must be reasonably small and 
was fixed at 0.0005tb. Initial configurations were gener- 
ated by putting the particles at random positions. For 
the rod systems, a parallel ncmatic initial configuration 
was adopted. After a long equilibration period of at least 
IOtb statistics were gathered (during a period of about 
15tb) and the time-dependent correlations were moni- 
tored during an interval of about 5tb- The latter turned 
out to be sufficiently large to reach the long-time limit. 



III. RESULTS FOR THE LONG-TIME 
SELF-DIFFUSION COEFFICIENTS 

Results for the long-time diffusion coefficient for 
spheres are shown in Fig. O A clear non-monotonic be- 
havior is observed. Both for very small and very high 
densities, the diffusivity comes very close to that of a sin- 
gle particle. At the highest density simulated (pa 3 = 25) 
the long-time diffusion constant has regained about 98 
% of its short-time value indicating that the diffusion 
has become virtually ideal in the high-density limit. The 
density pa 3 = 0.3 at which the diffusivity becomes mini- 
mal is in agreement with the results of Ref. [l6[ , and lies 
close to the density pa 3 = 0.23 for which the Gaussian 
core model displays its 'turning point' in the reentrant 
melting transition 

To compare our data with microscopic theory we have 
included the prediction from the Medina-Noyola theory 
for self-diffusion [l7j].The theory comprises an analysis 
of the effective Langevin equation of a tagged spherical 
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FIG. 3: Same as Fig. [2] for the Gaussian segment rods. 

colloid in a medium of interacting neighbor particles. In 
the absence of hydrodynamic interactions, the following 
expression for D\ is proposed: 

Dl=Dl(\+ P ~ Jdr[g 2 (r)-l] 2 ^) , (16) 

where the only input is the static pair correlation func- 
tion g 2 (r) which is obtained from the simulation. Al- 
though the theory is certainly not reliable from a quan- 
titative point of view, as we observe in Fig. [2j the non- 
monotonic behavior is clearly recovered. This suggests 
that there is a qualitative correspondence between the 
long-time diffusive behavior and the static correlations 
as embedded in g 2 (r). This becomes more explicit when 
we compare the diffusion data with the maximum ampli- 
tude in the pair correlation function, also shown in Fig.[5J 
The behavior for the rods is qualitatively the same as for 
spheres, see Fig. [3] Also here the translational diffu- 
sion constant varies non-monotonically with density and 
approaches the short-time limits at small and large den- 
sities. The only notable difference is that the fluid struc- 
ture is somewhat more pronounced here. As a result, the 
normalized diffusion constant reaches a minimum value 
that is smaller than that for spheres. Again, the the- 
ory of Medina-Noyola now taken with the center-of-mass 
pair correlations as an input overestimates the simulation 
data but shows the correct trend. 

Contrary to D\, the rotational counterpart in Fig. [4] 
seems to be weakly affected by the density. At the point 
of maximum fluid structure (per 3 = 0.12) the long-time 
rotational diffusivity has dropped to only about 90 % of 
the maximum i.e. short-time value. Moreover, an inves- 
tigation of W n (t) at minimum diffusivity shows that the 
mean-square orientational displacement does not depend 
on n. From this we may conclude that the rotational 
relaxation on the unit sphere is purely diffusive for long 
times. This behavior is not found in isotropic systems 
with unbounded rod potentials such as hard spherocylin- 
ders or Yukawa segment models [20|. The distinct dis- 




FIG. 4: (top) Mean-square displacement on the unit spheres 
W2(t) [see Eq. (|15[1 ] for Gaussian rods at various densities. 
The corresponding long-time rotational diffusion coefficients 
D — Dl/Dq are indicated where the number in brackets 
gives the error of the last digit, (bottom) W„(i) for n = 1, 2, 3 
for the system with density pa 3 = 0.12. 

crepancy between the density-dependence of the trans- 
lational and rotational diffusivity suggests that the cou- 
pling between orientational and translational degrees of 
freedom is very small for the systems of ultrasoft rods 
considered here. 

Some insight as to the status of the translation- rotation 
coupling can be gained by considering the effective inter- 
action of a rod pair in a spatially homogeneous fluid: 

v 2 (u>i,u>j) = J drijV2{rij ;u)t,u)j), (17) 

with Yij — Yj — r,-. Inserting Eq. and some algebra 
leads to «2 (<£>,-, &j) — const. This result holds for any 
bounded segment-segment potential and implies that, ir- 
respective of the rod aspect-ratio, all static orientational 
correlations are rendered zero by the random-phase ap- 
proximation for the excess free energy 0, [1[ . For a spa- 
tially uniform rod fluid the latter is simply proportional 
to a double orientational average of v 2 (fii,&j)- This will 
give the same outcome for any normalized orientational 
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distribution. Since the ideal free energy of a nematic fluid 
is always higher than that of the isotropic, the possibil- 
ity of a stable nematic state is fully excluded. Of course 
the above argument does not rule out a possible freezing 
transition occurring within an isotropic rod fluid. In fact, 
recent investigations for other soft rods such as parallel 
Gaussian-core particles [l5| and Yukawa rods [22j seem 
to point to a pronounced stability of columnar liquid- 
crystalline order in these systems. Finally, we remark 
that for soft rods with large aspect-ratios a phase tran- 
sition from an isotropic toward a nematic fluid may be 
possible at low densities where rod correlations are much 
better described by the Onsager functional [23| than the 
random-phase approximation. 

IV. CONCLUSIONS 

In conclusion, we have simulated the long-time self- 
diffusion in concentrated Brownian systems of rod-like 
and spherical particles which interact via a Gaussian core 
and are thus penetrable. As reflected by the statics, the 
system is getting ideal in the high-density limit where 
the random-phase approximation for the fluid structure 
becomes asymptotically exact. We think that the trends 
are independent of details in the interaction potentials 



provided that clustering [2| is avoided. 

We finish with a few remarks. First of all, one 
should consider the hydrodynamic interactions mediated 
by the solvent. These are neglected in the Brownian dy- 
namics simulations, but can be treated using more so- 
phisticated (and time-consuming) schemes like lattice- 
Boltzmann, Stokesian or the stochastic rotation dynam- 
ics [2J, |25|, l26|, l27|, [28|. Second, there is a need to de- 
rive microscopic models for the Brownian motion of stiff 
rods with soft interactions on the basis of the Smolu- 
chowski equation involving mode-coupling approxima- 
tions. These approaches then would go beyond a simple 
effective Medina-Noyola theory in treating explicitly the 
orientational degrees of freedom. Our simulation results 
may provide benchmark data to test these theories. As to 
the statics, it would be worthwhile to map out the freez- 
ing behavior of Gaussian segment rods in the regime of 
high density and low temperature. Finally, the existence 
of a stable nematic phase for Gaussian rods with large 
aspect-ratio remains an open question. 
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